A lo largo de esta práctica se ha realizado un caso práctico sobre el estudio de los datasets Red/White Wine Quality con el fin de aprender a identificar los datos relevantes de un proyecto analítico, así como utilizar las técnicas apropiadas para la integración, limpieza, validación y análisis de los datos.
A continuación, se muestran las competencias desarrolladas en esta práctica del Máster de Data Science.
* Capacidad de analizar un problema en el nivel de abstracción adecuado a cada situación y aplicar las habilidades y conocimientos adquiridos para abordarlo y resolverlo.
* Capacidad para aplicar las técnicas específicas de tratamiento de datos (integración, transformación, limpieza y validación) para su posterior análisis.
Como se especifica en el enunciado de esta práctica, los objetivos son los siguientes:
A menudo es inevitable, al producir una obra multimedia, hacer uso de recursos creados por terceras personas. Es por lo tanto comprensible hacerlo en el marco de una práctica de los estudios de Informática, Multimedia y Telecomunicación de la UOC, siempre y cuando esto se documente claramente y no suponga plagio en la práctica.
Por lo tanto, si se hicera uso de recursos ajenos, se presentará un documento detallando y especificando el nombre de cada recurso, su autor, el lugar dónde se obtuvo y su estatus legal: si la obra está protegida por el copyright o se acoge a alguna otra licencia de uso (Creative Commons, licencia GNU, GPL …). Como autores de la práctica nos aseguraremos que la licencia no impide específicamente su uso en el marco de la práctica. En caso de no encontrar la información correspondiente tendrá que asumir que la obra está protegida por copyright.
Como ya hemos comentado anteriormente, para el desarrollo de la práctica se ha analizado los datasets Red/White Wine Quality que podemos encontrar en UCI Machine Learning Repository. En estos ficheros cada registro representa un ensayo de calidad fisicoquímico sobre un determinado vino en base a una serie de atributos. Ambos ficheros presentan la misma configuración y el mismo número de atributos y registros como se muestra a continuación.
| Formato | CSV |
|---|---|
| Separador | ; |
| Header | True |
| Número de registros | 4898 |
| Número de atributos | 12 |
A continuación, se muestra una descripción de los atributos de los ficheros en estudio.
| Nombre | Descripción | Tipo |
|---|---|---|
| fixed acidity | Acidez fija | Numérico |
| volatile acidity | Acidez volátil | Numérico |
| citric acid | Ácido cítrico | Numérico |
| residual sugar | Azúcar residual | Numérico |
| chlorides | Cloruros | Numérico |
| free sulfur dioxide | Dióxido de azufre libre | Numérico |
| total sulfur dioxide | Dióxido de azufre total | Numérico |
| density | Densidad | Numérico |
| pH | PH | Numérico |
| sulphates | Sulfatos | Numérico |
| alcohol | Porcentaje de alcohol en volumen | Numérico |
| quality | Valor de calidad resultante del ensayo | Numérico |
El primer paso antes de realizar el análisis será unir ambos datasets para crear un único conjunto de datos.
El objetivo es realizar un análisis completo del dataset de modo que trabajemos todas las fases de la analítica avanzada.
Tipos de análisis
library(readr) # Providing a fast and friendly way to read rectangular data ('csv').
library(knitr) # A General-Purpose Package for Dynamic Report Generation in R.
library(tinytex) # Helper Functions to Install and Maintain 'TeX Live', and Compile 'LaTeX' Documents.
library(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntax. Build complex HTML or 'LaTeX' tables
library(tidyverse) # Collection of R packages designed for data science.
library(scales) # Graphical scales map data to aesthetics, and provide methods for automatically determining breaks and labels for axes
library(dplyr) # tool for working with data frame like objects, both in memory and out of memory.
library(tables) # Computes and displays complex tables of summary statistics
library(ggplot2) # A system for 'declaratively' creating graphics, based on "The Grammar of Graphics"
library(ggExtra) # Package for adding marginal histograms to ggplot2
library(magrittr) # Provides a mechanism for chaining commands with a new forward-pipe operator, %>%.
library(corrplot) # The corrplot package is a graphical display of a correlation matrix, confidence interval....
library(grid) # grid adds an nx by ny rectangular grid to an existing plot.
library(gridExtra) # Provides a number of user-level functions to work with "grid" graphics, notably to arrange multiple grid-based plots on a page, and draw tables.
library(visdat) # Create preliminary exploratory data visualisations of an entire dataset to identify problems or unexpected features using 'ggplot2'.
library(moments) # Compute skewness of a univariate distribution.
library(caret) # Pre-processing transformation (centering, scaling etc.) can be estimated from the training data and applied to any data set with the same variables.
library(nortest) # Performs the Lilliefors (Kolmogorov-Smirnov) test for the composite hypothesis of normality,
library(car) # VIF Regression: A Fast Regression Algorithm For Large Data
En primer lugar vamos a proceder a cargar ambos archivos. En este caso se han extraído directaemnte de la url, de modo que no tengan que pasar por nuestro dispositivo, optimizando el proceso. No obstante, se adjunta la línea de código comentada para cargar los archivos desde la carpeta raíz, por si hubiera alguna modificación de la url.
# Se carga el primer conjunto de datos
white_url <- "https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-white.csv"
whitewine <- read.csv(white_url, header = TRUE, sep = ";")
# whitewine <- read.csv(winequality-white.csv, header = TRUE, sep = ";")
# Se carga el segundo conjunto de datos
red_url <- "https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv"
redwine <- read.csv(white_url, header = TRUE, sep = ";")
#redwine <- read.csv(winequality-red.csv, header = TRUE, sep = ";")
#Añadimos una nueva columna con el color del vino antes de unir los dos conjuntos de datos.
whitewine$color<-"white"
redwine$color<-"red"
Unión de los conjuntos de datos
# Mediante la función rbind se unen (join) verticalmente los conjuntos de datos
wine <- rbind(redwine, whitewine)
Empezaremos haciendo un breve análisis de los datos ya que nos interesa tener una idea general de los datos que disponemos. Por ello, primero calcularemos las dimensiones de nuestra base de datos y analizaremos qué tipos de atributos tenemos.
Para empezar, calculamos las dimensiones mediante la función dim(). Obtenemos que disponemos de 9796 registros (filas) y 13 variables (columnas).
# Dimensiones de la base de datos
dim(wine)
## [1] 9796 13
# Estructura de los datos
str(wine)
## 'data.frame': 9796 obs. of 13 variables:
## $ fixed.acidity : num 7 6.3 8.1 7.2 7.2 8.1 6.2 7 6.3 8.1 ...
## $ volatile.acidity : num 0.27 0.3 0.28 0.23 0.23 0.28 0.32 0.27 0.3 0.22 ...
## $ citric.acid : num 0.36 0.34 0.4 0.32 0.32 0.4 0.16 0.36 0.34 0.43 ...
## $ residual.sugar : num 20.7 1.6 6.9 8.5 8.5 6.9 7 20.7 1.6 1.5 ...
## $ chlorides : num 0.045 0.049 0.05 0.058 0.058 0.05 0.045 0.045 0.049 0.044 ...
## $ free.sulfur.dioxide : num 45 14 30 47 47 30 30 45 14 28 ...
## $ total.sulfur.dioxide: num 170 132 97 186 186 97 136 170 132 129 ...
## $ density : num 1.001 0.994 0.995 0.996 0.996 ...
## $ pH : num 3 3.3 3.26 3.19 3.19 3.26 3.18 3 3.3 3.22 ...
## $ sulphates : num 0.45 0.49 0.44 0.4 0.4 0.44 0.47 0.45 0.49 0.45 ...
## $ alcohol : num 8.8 9.5 10.1 9.9 9.9 10.1 9.6 8.8 9.5 11 ...
## $ quality : int 6 6 6 6 6 6 6 6 6 6 ...
## $ color : chr "red" "red" "red" "red" ...
Una vez cargados los datos, vamos a verificar la transmisión de la información.
# Hacemos una primera visualización
head(wine)
## fixed.acidity volatile.acidity citric.acid residual.sugar chlorides
## 1 7.0 0.27 0.36 20.7 0.045
## 2 6.3 0.30 0.34 1.6 0.049
## 3 8.1 0.28 0.40 6.9 0.050
## 4 7.2 0.23 0.32 8.5 0.058
## 5 7.2 0.23 0.32 8.5 0.058
## 6 8.1 0.28 0.40 6.9 0.050
## free.sulfur.dioxide total.sulfur.dioxide density pH sulphates alcohol
## 1 45 170 1.0010 3.00 0.45 8.8
## 2 14 132 0.9940 3.30 0.49 9.5
## 3 30 97 0.9951 3.26 0.44 10.1
## 4 47 186 0.9956 3.19 0.40 9.9
## 5 47 186 0.9956 3.19 0.40 9.9
## 6 30 97 0.9951 3.26 0.44 10.1
## quality color
## 1 6 red
## 2 6 red
## 3 6 red
## 4 6 red
## 5 6 red
## 6 6 red
# Comprobamos el tipo de variable en los datos
sapply(wine,class)
## fixed.acidity volatile.acidity citric.acid
## "numeric" "numeric" "numeric"
## residual.sugar chlorides free.sulfur.dioxide
## "numeric" "numeric" "numeric"
## total.sulfur.dioxide density pH
## "numeric" "numeric" "numeric"
## sulphates alcohol quality
## "numeric" "numeric" "integer"
## color
## "character"
# Mostramos un resumen de los datos originales
summary(wine)
## fixed.acidity volatile.acidity citric.acid residual.sugar
## Min. : 3.800 Min. :0.0800 Min. :0.0000 Min. : 0.600
## 1st Qu.: 6.300 1st Qu.:0.2100 1st Qu.:0.2700 1st Qu.: 1.700
## Median : 6.800 Median :0.2600 Median :0.3200 Median : 5.200
## Mean : 6.855 Mean :0.2782 Mean :0.3342 Mean : 6.391
## 3rd Qu.: 7.300 3rd Qu.:0.3200 3rd Qu.:0.3900 3rd Qu.: 9.900
## Max. :14.200 Max. :1.1000 Max. :1.6600 Max. :65.800
## chlorides free.sulfur.dioxide total.sulfur.dioxide
## Min. :0.00900 Min. : 2.00 Min. : 9.0
## 1st Qu.:0.03600 1st Qu.: 23.00 1st Qu.:108.0
## Median :0.04300 Median : 34.00 Median :134.0
## Mean :0.04577 Mean : 35.31 Mean :138.4
## 3rd Qu.:0.05000 3rd Qu.: 46.00 3rd Qu.:167.0
## Max. :0.34600 Max. :289.00 Max. :440.0
## density pH sulphates alcohol
## Min. :0.9871 Min. :2.720 Min. :0.2200 Min. : 8.00
## 1st Qu.:0.9917 1st Qu.:3.090 1st Qu.:0.4100 1st Qu.: 9.50
## Median :0.9937 Median :3.180 Median :0.4700 Median :10.40
## Mean :0.9940 Mean :3.188 Mean :0.4898 Mean :10.51
## 3rd Qu.:0.9961 3rd Qu.:3.280 3rd Qu.:0.5500 3rd Qu.:11.40
## Max. :1.0390 Max. :3.820 Max. :1.0800 Max. :14.20
## quality color
## Min. :3.000 Length:9796
## 1st Qu.:5.000 Class :character
## Median :6.000 Mode :character
## Mean :5.878
## 3rd Qu.:6.000
## Max. :9.000
Factorizamos la variable “color” para futuros análisis.
wine$color<-as.factor(wine$color)
Se eliminan los registros duplicados si los hubiese
wine <- wine[!duplicated(wine), ]
dim(wine)
## [1] 7922 13
Comprobamos sie existen valores NA
vis_miss(wine)
Todos los valores están presentes. No existen valores NA.
Como podemos ver, todas las variables han sido identificadas con el tipo correspondiente y no se aprecian anomalías a simple vista en los datos.
attach(wine)
#Iniciamos el cuadrante para el gráfico
par(mfrow=c(1,2))
# Representamos histogramas y boxplot para las variables numéricas
hist(fixed.acidity , main="fixed.acidity histogram")
boxplot(fixed.acidity, main="fixed.acidity boxplot")
hist(volatile.acidity , main="volatile.acidity histogram")
boxplot(volatile.acidity , main="volatile.acidity boxplot")
hist(citric.acid , main="citric.acid histogram")
boxplot(citric.acid , main="citric.acid boxplot")
hist(residual.sugar , main="residual.sugar histogram")
boxplot(residual.sugar , main="residual.sugar boxplot")
hist(chlorides , main="chlorides histogram")
boxplot(chlorides , main="chlorides boxplot")
hist(free.sulfur.dioxide , main="free.sulfur.dioxide histogram")
boxplot(free.sulfur.dioxide , main="free.sulfur.dioxide boxplot")
hist(total.sulfur.dioxide , main="total.sulfur.dioxide histogram")
boxplot(total.sulfur.dioxide , main="total.sulfur.dioxide boxplot")
hist(density , main="density histogram")
boxplot(density , main="density boxplot")
hist(sulphates , main="sulphates histogram")
boxplot(sulphates , main="sulphates boxplot")
hist(alcohol , main="alcohol histogram")
boxplot(alcohol , main="alcohol boxplot")
hist(quality , main="quality histogram")
boxplot(quality , main="quality boxplot")
# Representamos gráfico de barras de las variables nominales
barplot(table(color), main="Color histogram")
detach(wine)
Todas las variables presentan outliers.
Comprobamos a continuación la asimetría de las variables verificar si los datos se distribuyen normalmente o no.
attach(wine)
skewness(fixed.acidity)
## [1] 0.6958366
skewness(volatile.acidity)
## [1] 1.640459
skewness(citric.acid)
## [1] 1.310105
skewness(residual.sugar)
## [1] 1.333134
skewness(chlorides)
## [1] 4.967194
skewness(free.sulfur.dioxide)
## [1] 1.566087
skewness(total.sulfur.dioxide)
## [1] 0.4566267
skewness(density)
## [1] 1.272836
skewness(pH)
## [1] 0.4552843
skewness(sulphates)
## [1] 0.9374981
skewness(alcohol)
## [1] 0.4505259
skewness(quality)
## [1] 0.1119616
detach(wine)
Regla de la asimetría:
Fuente:*GoodData.Documentation
Se aplica transformación de Boxcox y volvemos a comprobar la asimetría.
Fuente: *Box-Cox Transformations
wine_procesado <- preProcess(wine[,1:11], c("BoxCox", "center", "scale"))
new_wine <- data.frame(trans = predict(wine_procesado, wine))
colnames(new_wine)
## [1] "trans.fixed.acidity" "trans.volatile.acidity"
## [3] "trans.citric.acid" "trans.residual.sugar"
## [5] "trans.chlorides" "trans.free.sulfur.dioxide"
## [7] "trans.total.sulfur.dioxide" "trans.density"
## [9] "trans.pH" "trans.sulphates"
## [11] "trans.alcohol" "trans.quality"
## [13] "trans.color"
attach(new_wine)
skewness(trans.fixed.acidity)
## [1] 0.1131816
skewness(trans.volatile.acidity)
## [1] 0.1719291
skewness(trans.citric.acid)
## [1] 1.310105
skewness(trans.residual.sugar)
## [1] -0.05828706
skewness(trans.chlorides)
## [1] -0.1496725
skewness(trans.free.sulfur.dioxide)
## [1] 0.0950706
skewness(trans.total.sulfur.dioxide)
## [1] 0.008934387
skewness(trans.density)
## [1] 1.100076
skewness(trans.quality)
## [1] 0.1119616
skewness(trans.pH)
## [1] 0.0002882357
skewness(trans.sulphates)
## [1] -0.01873998
skewness(trans.alcohol)
## [1] 0.03930492
detach(new_wine)
Ahora la mayoría las distribuciones son aproximadamente simétricas (la asimetría está entre -0,5 y 0,5). Las únicas variables que muestran aún una asimestría son:
* trans.citric.acid * trans.density
Se va a realizar ejemplarmente una comprobación de la normalidad para la variable trans.citric.acid.
# Para comprobar la normalidad de la variable trans.citric.acid se va a usar:
# 1. Gráfico de cuantiles teóricos (Gráficos Q-Q): Consiste en comparar los cuantiles de la distribución observada con los cuantiles teóricos de una distribución normal con la misma media y desviación estándar que los datos. Cuanto más se aproximen los datos a una normal, más alineados están los puntos entorno a la recta.
# 2. Un histograma
# 3. Diagrama de cajas
par(mfrow = c(1, 3))
hist(new_wine$trans.citric.acid, las=1, main="Distribución normal", font.main=1, ylab="Frecuencia", xlab = "trans.citric.acid")
qqnorm(new_wine$trans.citric.acid, las=1, pch=18, main="Simetria", font.main=1, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales")
qqline(new_wine$trans.citric.acid)
boxplot(new_wine$trans.citric.acid, las=1, main="Valores extremos", font.main=1, xlab="Datos de distribution normal", horizontal=F)
No se trata de una distribución normal ya que presenta una asimetría en ambas colas, aunque la cola a la derecha de la media es más larga que la izquierda.
Para asegurarnos completamente se realiza la pruena de Lilliefors (Kolmogorov-Smirnov) para muestras grandes (n>50).
El test Lilliefors asume que la media y varianza son desconocidas, estando especialmente desarrollado para contrastar la normalidad.
Es la alternativa al test de Shapiro-Wilk cuando el número de observaciones es mayor de 50. La función lillie.test() del paquete nortest permite aplicarlo.
Planteamos la Hipótesis:
\(\ H_0:\) La muestra proviene de una distribución normal.
\(\ H_1:\) La muestra no proviene de una distribución normal.
Nivel de significancia: El nivel de significancia que se trabajará es de 0.05. \(\alpha=0.05\)
Criterio de decisón:
Si \(p<\alpha\) se rechaza \(\ H_0:\)
Si \(p>\alpha\) no se rechaza \(\ H_0:\)
lillie.test(x=new_wine$trans.citric.acid)
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: new_wine$trans.citric.acid
## D = 0.10992, p-value < 2.2e-16
Dado que \(p-value<\alpha\) se rechaza \(\ H_0:\). La muestra no proviene de una distribución normal. Pero, puesto que el tamaño de la muestra es mayor a 30, según el teorema del límite central, asumimos la normalidad.
Procedemos a eliminar los valores outliers encontardos en los gráficos boxplot anteriores.
# Mostramos los valores identificados como outliers
attach(new_wine)
f_a_out<-boxplot.stats(trans.fixed.acidity)$out
v_a_out<-boxplot.stats(trans.volatile.acidity)$out
c_a_out<-boxplot.stats(trans.citric.acid)$out
r_s_out<-boxplot.stats(trans.residual.sugar)$out
ch_out<-boxplot.stats(trans.chlorides)$out
f_s_d_out<-boxplot.stats(trans.free.sulfur.dioxide)$out
t_s_d_out<-boxplot.stats(trans.total.sulfur.dioxide)$out
d_t_out<-boxplot.stats(trans.density)$out
ph_out<-boxplot.stats(trans.pH )$out
sh_out<-boxplot.stats(trans.sulphates)$out
ac_out<-boxplot.stats(trans.alcohol)$out
qa_out<-boxplot.stats(trans.quality)$out
detach(new_wine)
Las variables trans.alcohol, trans.residual.sugar y trans.color no presentan outliers.
# Eliminamos las instancias con outlier
new_wine<-new_wine[-which(new_wine$trans.fixed.acidity %in% f_a_out),]
new_wine<-new_wine[-which(new_wine$trans.volatile.acidity %in% v_a_out),]
new_wine<-new_wine[-which(new_wine$trans.citric.acid %in% c_a_out),]
new_wine<-new_wine[-which(new_wine$trans.chlorides %in% ch_out),]
new_wine<-new_wine[-which(new_wine$trans.free.sulfur.dioxide %in% f_s_d_out),]
new_wine<-new_wine[-which(new_wine$trans.total.sulfur.dioxide %in% t_s_d_out),]
new_wine<-new_wine[-which(new_wine$trans.density %in% d_t_out),]
new_wine<-new_wine[-which(new_wine$trans.pH %in% ph_out),]
new_wine<-new_wine[-which(new_wine$trans.sulphates %in% sh_out),]
new_wine<-new_wine[-which(new_wine$trans.quality %in% qa_out),]
Tras haber eliminado los valores atípicos, revisamos la correlación entre las distintas variables. El objetivo es centrar el análisis descriptivo en aquellas variables que presentan correlación.
#Se omite la variable color para pocer representar la matriz de correlación
corrplot(cor(new_wine[1:12]))
Las correlaciones positivas se muestran en azul y las negativas en rojo. La intensidad del color y el tamaño del círculo son proporcionales a los coeficientes de correlación. En la parte derecha del correlograma, el color de la leyenda muestra los coeficientes de correlación y los colores correspondientes.
Observando los resultados se verifica que:
1. La calidad del vino no se encuentra correlacionada con las variables “citric.acid”, “free.sulfur.dioxide”, los “sulphates” y “color”. 2. Es coherente la alta correlación positiva entre free.sulfur.dioxide y total.sulfur.dioxide. A mayor cantidad de sulfuro más sulfuro libre hay. 3. Es coherente la alta correlación negativa entre fixed.acidity y PH. Cuanto mayor es el PH en la escala de ácido-base más básico es el compuesto.Por lo tanto a valores altos de PH menor es la acidez. 4. Es coherente la alta correlación negativa entre el azucar y el alcohol.El alcohol es el resultado de la oxidación del azucar, por tanto si hay menos azucar el grado de alcohol aumenta.
A continuación se estudian visaulamente las posibles correlaciones entre:
# Calculamos la correlación y la mostramos por pantalla
RS_D_cor = cor(new_wine$trans.residual.sugar,new_wine$trans.density)
print(RS_D_cor)
## [1] 0.76787
# Generamos el scatter plot y añadimos la linea de regresión para una mejor visualización
plot(new_wine$trans.residual.sugar,new_wine$trans.density, main="residual.sugar vs density",xlab="residual.sugar", ylab="density", pch=19)
abline(lm(new_wine$trans.density~new_wine$trans.residual.sugar), col="red") # regression line (y~x)
Como podemos apreciar en la gráfica existe una fuerte relación lineal entre ambas variables y todos los puntos se ajustan bastante bien a la linea. El resultado de la correlación, 0.75, refuerza lo comentado anteriormente.
# Calculamos la correlación y la mostramos por pantalla
D_A_cor = cor(new_wine$trans.density,new_wine$trans.alcohol)
print(D_A_cor)
## [1] -0.806588
# Generamos el scatter plot y añadimos la linea de regresión para una mejor visualización
plot(new_wine$trans.density,new_wine$trans.alcohol, main="alcochol vs density",xlab="alcochol", ylab="density", pch=19)
abline(lm(new_wine$trans.alcohol~new_wine$trans.density), col="red") # regression line (y~x)
Al igual que antes, existe una fuerte relación lineal entre ambas variables y los puntos se ajustan a la linea. El resultado de la correlación, -0.80, refuerza lo comentado anteriormente.
# Calculamos la correlación y la mostramos por pantalla
par(mfrow = c(1,2))
for (i in c(1:11)) {
plot(new_wine[, i], jitter(new_wine[, "trans.quality"]), xlab = names(new_wine)[i],
ylab = "trans.quality", col = "firebrick", cex = 0.8, cex.lab = 1.3)
abline(lm(new_wine[, "trans.quality"] ~ new_wine[ ,i]), lty = 2, lwd = 2)
}
par(mfrow = c(1, 1))
Las variables trans.alcohol, trans.sulphates y trans.pH parecen tener una correlación positiva con trans.quality, mientras que ltrans.volatile.acidity, trans.chlorides y trans.total.sulfur.dioxidel parecen tener una relación negativa con la calidad.
Dado que nuestro dataset contiene los valores para los vinos blancos y tintos, vamos a realizar un contraste de hipótesis para la diferencia de medias (trans.citric.acid) de las dos muestras. Planteamos la sigueinte hipótesis:
¿Podemos aceptar con un nivel de confianza del 95% que los vinos blancos (color) tienen una acidez (trans.citric.acid) que supera a la acidez de los vinos tintos (color)?.
Planteamos la Hipótesis nula y alternativa \[ \left\{ \begin{array}{ll} H_{0}: & \mu_{si}=\mu_{no}\\ H_{1}: & \mu_{si}>\mu_{no} \end{array} \right.\\ \ Hipótesis\ unilateral\ a \ la \ dereccha \]
EL método que se aplica es un contraste de hipótesis para dos muestras independientes. El test es paramétrico, asumiendo el teorema del límite central. Es un test unilateral, con varianzas desconocidas, pero iguales, ya que el test de varianzas no rechaza la hipótesis nula de varianzas iguales.
#Test de igualdad de varianzas
#H0: F (ratio de varianzas) = 1
#H1: F diferente de 1
var.test(new_wine$trans.citric.acid[new_wine$trans.color=="white"], new_wine$trans.citric.acid[new_wine$trans.color=="red"],
alternative = "greater")
##
## F test to compare two variances
##
## data: new_wine$trans.citric.acid[new_wine$trans.color == "white"] and new_wine$trans.citric.acid[new_wine$trans.color == "red"]
## F = 1, num df = 3242, denom df = 3242, p-value = 0.5
## alternative hypothesis: true ratio of variances is greater than 1
## 95 percent confidence interval:
## 0.9438529 Inf
## sample estimates:
## ratio of variances
## 1
A continuación se lleva a cabo el contraste de hipótesis:
#Se define una función myttest.
#x1, x2: muestras
#CL: confidence level (nivel de confianza)
#equalvar TRUE/FALSE
#alternative: "bilateral", "less" (x1 less than x2), "greater" (x1 greater than x2).
myttest <- function( x1, x2, CL=0.95,equalvar=TRUE, alternative="bilateral" ){ # z test
mean1<-mean(x1)
n1<-length(x1)
sd1<-sd(x1)
mean2<-mean(x2)
n2<-length(x2)
sd2<-sd(x2)
if (equalvar==TRUE){
s <-sqrt( ((n1-1)*sd1^2 + (n2-1)*sd2^2 )/(n1+n2-2) )
Sb <- s*sqrt(1/n1 + 1/n2)
df<-n1+n2-2
}
else{ #equalvar==FALSE
Sb <- sqrt( sd1^2/n1 + sd2^2/n2 )
denom <- ( (sd1^2/n1)^2/(n1-1) + (sd2^2/n2)^2/(n2-2))
df <- ((sd1^2/n1 + sd2^2/n2)^2) / denom
}
alfa <- (1-CL)
t<- (mean1-mean2) / Sb
if (alternative=="bilateral"){
tcritical <- qt( alfa/2, df, lower.tail=FALSE ) #two sided
pvalue<-pt( abs(t), df, lower.tail=FALSE )*2 #two sided
}
else if (alternative=="less"){
tcritical <- qt( alfa, df, lower.tail=TRUE )
pvalue<-pt( t, df, lower.tail=TRUE )
}
else{ #(alternative=="greater")
tcritical <- qt( alfa, df, lower.tail=FALSE )
pvalue<-pt( t, df, lower.tail=FALSE )
}
#Guardamos en resultado en un data frame
info<-data.frame(t,df,tcritical,pvalue)
info %>% kable() %>% kable_styling()
return (info)
}
Se visualizan tanto el estadístico como el valor p:
info<-myttest( new_wine$trans.citric.acid[new_wine$trans.color=="white"],
new_wine$trans.citric.acid[new_wine$trans.color=="red"],
alternative="less",
equalvar="TRUE")
info
## t df tcritical pvalue
## 1 0 6484 -1.645089 0.5
Dado que hemos aplicado una función propia comprobamos el resultado con la fución original de R.
#Comparación con t.test de R
t.test( new_wine$trans.citric.acid[new_wine$trans.color=="white"], new_wine$trans.citric.acid[new_wine$trans.color=="red"], alternative = "less",
var.equal = TRUE)
##
## Two Sample t-test
##
## data: new_wine$trans.citric.acid[new_wine$trans.color == "white"] and new_wine$trans.citric.acid[new_wine$trans.color == "red"]
## t = 0, df = 6484, p-value = 0.5
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
## -Inf 0.02953236
## sample estimates:
## mean of x mean of y
## -0.07990945 -0.07990945
El valor p es 0.5, mayor que el nivel de significación es 0.05 (95%). No podemos rechazar la hipótesis nula de igualdad de acidez media entre los vinos blancos y tintos. En otras palabras, no podemos asumir que la media de acidez de vinos blancos sea mayor a la de vinos tintos.
Realizad los cálculos manualmente del test. Para ello, se recomienda construir una función my.chisq que realice estos cálculos y que podáis reusar más adelante.
# Realizamos el test con la función chisq.test para comprobar los resultados.
chisq.test(new_wine$trans.color,new_wine$trans.quality)
##
## Pearson's Chi-squared test
##
## data: new_wine$trans.color and new_wine$trans.quality
## X-squared = 0, df = 3, p-value = 1
El valor crítico es mayor que el chi-cuadrado, por tanto, no podemos rechazar la hipótesis nula de que las variables son independientes.
Del mismo modo, el p valor es mayor al nivel de confianza 0.05 por tanto también debemos aceptar la hipótesis nula y podemos afirmar que no hay una dependencia entre el color y si la calidad del vino con un 95% de confianza.
Vamos a aplicar un modelo de regresión lineal multivariable. El objetivo consistirá en predecir la calidad del vino en función de las otras variables.
Antes de elegir el mejor modelo lineal vamos a dividir el conjunto de datos en dos:
1. Set de entrenamiento
2. Set de prueba
set.seed(100)
#new_wine$trans.color<-as.numeric(new_wine$trans.color)
trainingRowIndex <- sample(1:nrow(new_wine), 0.8*nrow(new_wine))
new_wine_train <- new_wine[trainingRowIndex, ]
new_wine_test <- new_wine[-trainingRowIndex, ]
Realizamos el primer modelo
# Incluimos todas las variables
modelo_1 <- lm(trans.quality ~ . , new_wine_train)
summary(modelo_1)
##
## Call:
## lm(formula = trans.quality ~ ., data = new_wine_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.44993 -0.43246 -0.00713 0.45640 2.18459
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.812321 0.012632 460.132 < 2e-16 ***
## trans.fixed.acidity 0.058439 0.014282 4.092 4.35e-05 ***
## trans.volatile.acidity -0.143114 0.010565 -13.545 < 2e-16 ***
## trans.citric.acid 0.035017 0.013096 2.674 0.007523 **
## trans.residual.sugar 0.261757 0.022333 11.720 < 2e-16 ***
## trans.chlorides -0.077566 0.014166 -5.476 4.57e-08 ***
## trans.free.sulfur.dioxide 0.143246 0.012859 11.140 < 2e-16 ***
## trans.total.sulfur.dioxide -0.049529 0.014211 -3.485 0.000496 ***
## trans.density -0.342990 0.037590 -9.125 < 2e-16 ***
## trans.pH 0.095568 0.012167 7.855 4.84e-15 ***
## trans.sulphates 0.073967 0.010040 7.368 2.01e-13 ***
## trans.alcohol 0.164079 0.023736 6.913 5.33e-12 ***
## trans.colorwhite 0.005206 0.017739 0.293 0.769169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6388 on 5175 degrees of freedom
## Multiple R-squared: 0.2883, Adjusted R-squared: 0.2866
## F-statistic: 174.7 on 12 and 5175 DF, p-value: < 2.2e-16
Todas las variables presentan significación muy alta, esto puede ser consecuencia de multicolinealidad entre las variables. Para detectarla y minimizarla vamos a seguir los siguientes pasos:
Multicolinealidad: Se dice que existe multicolinealidad entre las variables explicativas cuando existe algún tipo de dependencia lineal entre ellas, o lo que es lo mismo, si existe una fuerte correlación entre las mismas. La correlación no solamente se refiere a las distintas variables dos a dos, sino a cualquier de ellas con cualquier grupo de las restantes. Por esta razón no es suficiente (aunque sí necesaria) que en la matriz de correlaciones bivariadas haya correlaciones altas.
El principal inconveniente de la multicolinealidad consiste en que se incrementan la varianza de los coeficientes de regresión estimados hasta el punto que resulta prácticamente imposible establecer su significación estadística, ya que como se sabe, el valor de t para un determinado coeficiente de regresión es el valor de dicho coeficiente dividido por su desviación tipica. Si este es grande, el valor de t será bajo y no llegara a la significación.
Hay varios procedimientos para detectar la multicolinealidad entre los predictores. El primero de ellos, basado en la correlación múltiple de un determinado regresor con los restantes se denomina Tolerancia de dicho regresor. Su valor es \(\ 1-R^2_i\). Siendo \(\ R^2_i\) la correlación múltiple al cuadrado de dicho regresor con los restantes.
Para que haya multicolinealidad dicha correlación ha de ser alta, o lo que es lo mismo la tolerancia baja.
Otro índice relacionado con éste y que nos da una idea del grado de aumento de la varianza se denomina Factor de Inflación de la Varianza, y es precisamente el recíproco de la tolerancia. Su valor es: \(\ VIF=1/(1-R^2_y)\). Para que no haya multicolinealidad el denominador tiene que valer cerca de la unidad, por tanto un poco más de 1 el valor de VIF. Cuanto mayor sea de este valor mayor multicolinealidad habrá y por tanto más inestable es el modelo lineal; en otras palabras, es conveniente eliminar del modelo las variables con un factor VIF grande.
Si todos los VIF son 1, no hay multicolinealidad, pero si algunos VIF son mayores que 1, los predictores están correlacionados. Cuando un VIF es > 5, el coeficiente de regresión para ese término no se estima adecuadamente.
Fuentes: REGRESIÓN LINEAL MÚLTIPLE
Multicolinealidad
Analizamos la multicolinealidad mediante el análisis del valor VIF
vif(modelo_1)
## trans.fixed.acidity trans.volatile.acidity
## 2.057763 1.178893
## trans.citric.acid trans.residual.sugar
## 1.133096 6.259537
## trans.chlorides trans.free.sulfur.dioxide
## 1.462590 1.848040
## trans.total.sulfur.dioxide trans.density
## 2.372484 16.362994
## trans.pH trans.sulphates
## 1.652030 1.108447
## trans.alcohol trans.color
## 6.666656 1.000248
Se observa multicolinealidad en la variable trans.density y se procede a eliminarla del modelo.
Ajustamos el modelo de nuevo
# excluímos la variable trans.density
modelo_2 <- lm(trans.quality ~ . -trans.density, new_wine_train)
summary(modelo_2)
##
## Call:
## lm(formula = trans.quality ~ . - trans.density, data = new_wine_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.41233 -0.42608 -0.00765 0.44988 2.21517
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.814549 0.012729 456.779 < 2e-16 ***
## trans.fixed.acidity -0.021857 0.011338 -1.928 0.053941 .
## trans.volatile.acidity -0.137441 0.010631 -12.929 < 2e-16 ***
## trans.citric.acid 0.024103 0.013145 1.834 0.066763 .
## trans.residual.sugar 0.082779 0.010763 7.691 1.74e-14 ***
## trans.chlorides -0.094106 0.014161 -6.646 3.33e-11 ***
## trans.free.sulfur.dioxide 0.151395 0.012929 11.709 < 2e-16 ***
## trans.total.sulfur.dioxide -0.061560 0.014262 -4.316 1.62e-05 ***
## trans.pH 0.040241 0.010632 3.785 0.000156 ***
## trans.sulphates 0.055080 0.009902 5.563 2.79e-08 ***
## trans.alcohol 0.349413 0.012378 28.228 < 2e-16 ***
## trans.colorwhite 0.004975 0.017880 0.278 0.780843
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6438 on 5176 degrees of freedom
## Multiple R-squared: 0.2768, Adjusted R-squared: 0.2753
## F-statistic: 180.1 on 11 and 5176 DF, p-value: < 2.2e-16
Analizamos la multicolinealidad mediante el análisis del valor VIF
vif(modelo_2)
## trans.fixed.acidity trans.volatile.acidity
## 1.276563 1.174811
## trans.citric.acid trans.residual.sugar
## 1.123644 1.431112
## trans.chlorides trans.free.sulfur.dioxide
## 1.438643 1.839125
## trans.total.sulfur.dioxide trans.pH
## 2.352063 1.241750
## trans.sulphates trans.alcohol
## 1.061328 1.784764
## trans.color
## 1.000246
Los valores VIF son menores de 5 para todas las variables, en otras palabras la colinealidad es muy baja con el modelo 2.
Ajustamos un modelo nuevo con las variables que presentan significación:
modelo_3 <- lm(trans.quality ~ trans.volatile.acidity +trans.free.sulfur.dioxide +trans.alcohol , new_wine_train)
summary(modelo_3)
##
## Call:
## lm(formula = trans.quality ~ trans.volatile.acidity + trans.free.sulfur.dioxide +
## trans.alcohol, data = new_wine_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.43180 -0.43398 -0.01586 0.47378 2.08740
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.816748 0.009104 638.89 <2e-16 ***
## trans.volatile.acidity -0.140007 0.009994 -14.01 <2e-16 ***
## trans.free.sulfur.dioxide 0.145472 0.010077 14.44 <2e-16 ***
## trans.alcohol 0.382024 0.009762 39.13 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6538 on 5184 degrees of freedom
## Multiple R-squared: 0.2532, Adjusted R-squared: 0.2527
## F-statistic: 585.8 on 3 and 5184 DF, p-value: < 2.2e-16
Los valores VIF son menores de 1.1 para todas las variables, en otras palabras la no hay colinealidad en el modelo 3.
vif(modelo_3)
## trans.volatile.acidity trans.free.sulfur.dioxide
## 1.006974 1.083543
## trans.alcohol
## 1.076535
Predecimos valores utilizando el set de prueba
# Se redondea el resultado para eliminar los decimales en la predicción
modelo_3_Pred = round(predict(modelo_3, newdata=new_wine_test))
head(modelo_3_Pred)
## 12 16 25 34 58 61
## 5 6 5 6 5 5
Creamos la matriz de confusión para comprobar la precisión del modelo. Se enfrentan los valores predichos por el modelo con los valores reales del set de prueba.
confusionMatrix(table(modelo_3_Pred, new_wine_test$trans.quality))
## Confusion Matrix and Statistics
##
##
## modelo_3_Pred 4 5 6 7
## 4 2 0 0 0
## 5 21 159 100 10
## 6 19 216 502 243
## 7 0 0 11 15
##
## Overall Statistics
##
## Accuracy : 0.5223
## 95% CI : (0.4948, 0.5498)
## No Information Rate : 0.4723
## P-Value [Acc > NIR] : 0.00017
##
## Kappa : 0.1689
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 4 Class: 5 Class: 6 Class: 7
## Sensitivity 0.047619 0.4240 0.8189 0.05597
## Specificity 1.000000 0.8581 0.3022 0.98932
## Pos Pred Value 1.000000 0.5483 0.5122 0.57692
## Neg Pred Value 0.969136 0.7857 0.6509 0.80110
## Prevalence 0.032357 0.2889 0.4723 0.20647
## Detection Rate 0.001541 0.1225 0.3867 0.01156
## Detection Prevalence 0.001541 0.2234 0.7550 0.02003
## Balanced Accuracy 0.523810 0.6410 0.5606 0.52265
Nuestro modelo tiene una precisión del 52.23%
Diagnosticamos el modelo a continuación
#Revisamos los residuos para asegurar que están distribuidos normalmente
hist(modelo_3$residuals)
#Comprobamos la media de los residuos para asegurarnos que ésta se encuentra cerca de 0
mean(modelo_3$residuals)
## [1] 1.041235e-17
#Visualizamos los diagramas
layout(matrix(c(1,2,3,4),2,2))
plot(modelo_3)
Tras examinar los resultados podemos afirmar: * Se satisface la condición de linealidad.
* El análisis gráfico confirman la normalidad.
* No hay evidencias de falta de homocedasticidad.
* No hay multicolinialidad.
Para finalizar se genera un archivo con el conjunto de datos procesado para análisis posteriores.
# Guardamos los datos
write.csv(new_wine, file="Wine_clean.csv",row.names=F)
A continuación se muestra la dirección de acceso al repositorio, donde podemos encontrar el código fuente, el archivo original, el archivo preprocesado y el documento html (RMarkdown).
Tanto la la tabla de contribuciones como las referncias utilizadass durante esta práctica, se encuentran disponibles en la wiki del repositorio.